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The reheating dynamics after the inflation induced by i?^-corrected f{R) model is considered. To 
avoid the complexity of solving the fourth order equations, we analyze the inflationary and reheating 
dynamics in the Einstein frame and its analytical solutions are derived. We also perform numerical 
calculation including the backreaction from the particle creation and compare the results with the 
analytical solutions. Based on the results, observational constraints on the model are discussed. 
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' Seeking the physical origin of two accelerated expansion regimes of the Universe, namely, the primordial inflation 
and the present cosmic acceleration, is one of the most important theoretical challenges of cosmology today. These 
unknown physical origins are referred by the primordial dark energy (DE) and the present DE. Various theoretical 
_ . models have been proposed to accelerate the cosmological expansion, 
'w , Among those theoretical models, f{R) gravity is a simple and nontrivial generalization of General Relativity. For 
(~| ' a recent review, see Refs. [1, 2]. It introduces a function f{R) in the action, where R is Ricci curvature. This 
Qh! additional degree of freedom plays a role of a scalar field, which is called scalaron, and it can cause the accelerated 
Q ' expansion of the Universe. The original idea recognized as R^ inflation model was proposed in Ref. [3], where de 
\ Sitter expansion was derived as a solution for the Einstein equation with quantum one-loop correction. After the 

■ accelerating expansion, the particles are gravitationally created through the oscillation of the scalaron, and it leads 
' radiation dominated Universe [4-8]. The R^ model predicts an almost scale-invariant spectrum, whose scalar and 

tensor components are consistent with recent observational data (Sec e.g., [9]). 
' Later, the R^ gravity was extended to a general function of R, namely f{R), to describe the late time cosmic 
^ ' acceleration. After some early challenges, the viable f{R) models were proposed [10-12], which can realize stable 

■ matter-dominated regime and subsequent late time acceleration. In these models, the expansion history of the Universe 
I is close to that in the A-cold dark matter (ACDM) model. From the model selection point of view, the key to distinguish 
. the models is focusing on small deviation from the ACDM model. As for background quantity, the equation of state 

" I ' parameter wde well characterizes the models. While it remains constant wde = — 1 in the ACDM model, it is time 
■rj" dependent in the f{R) gravity and even crosses the phantom divide at redshift z ^ 0{1) [10, 13, 14]. On the other 
' hand, the growth of the matter density fluctuations is also useful to measure the deviation. Since in f{R) gravity the 
, effective gravitational constant depends on time and scale, the matter power spectrum is enhanced [10, 12, 15-20]. 
This enhancement not only measure the deviation but also provides another interesting consequence: f{R) gravity 
admits larger neutrino mass. This is because massive neutrinos suppress evolution of matter fluctuations by free 
• streaming, which cancels the enhancement in f{R) gravity. As a result, f{R) gravity allows larger neutrino mass up 
I to ~ 0.5 eV [21]. The constraint on sterile neutrino mass is also relaxed up to ^ 1 eV, which is consistent with recent 
d . experiments [22]. Other distinguishable features of f{R) gravity would be imprinted on cosmological gravitational 
waves [23-25] . Future pulsar timing experiments and gravitational- wave detectors will be able to probe them directly 
and test gravity theories [26-29] . 

However, the above viable f{R) models still have theoretical problems [12]. If we start from some initial condition 
and calculate back to the past, the scalaron mass diverges quickly and the scalaron oscillates rapidly [30]. Another 
problem is that the Ricci scalar also diverges at the past even if we include nonlinear effect [31]. This curvature 
singularity was also pointed out in Refs. [32, 33]. To solve these problems, i?^-corrected f{R) model has been 
proposed [34]. This model is constructed from the late time acceleration part and R^ term. Consequently, the 
reheating followed by inflation in this model is different from that in R'^ model, i.e., the scalaron does not oscillate 
harmonically. Of course, we can construct the other specific functional forms that avoid singularities and describe 
both the primordial DE and the present DE. However, these functions belong to the same class and their behavior are 
similar [35]. In this sense, it is worth to study one specific model in detail as an example of such a class of extended 
f{R) models. 

The aim of this paper is to investigate the evolution during the inflation and reheating regimes of the i?^-corrected 
f{R) model in detail. In the previous work [34], the reheating dynamics has been analyzed in the Jordan frame. 
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However, the field equation is fourth order differential equation, which makes the physical interpretation unclear. 
In addition, in their analysis, the inflation and the reheating are separately solved with different initial conditions. 
Therefore, we focus on the evolution of the scalaron rolling on a potential in the Einstein frame, which clarifies the 
physical picture and allows us to understand the dynamics intuitively. We start from a certain initial condition 
imposed during inflationary regime and numerically solve the transition from the inflation to the reheating and the 
following reheating dynamics. Thus, our analysis is more accurate than the previous work [34] and is complementary 
to the analysis in the Jordan frame. 

This paper is organized as follows. In Sec. II, we review the basic equations in f{R) gravity in both the Jordan 
frame and the Einstein frame. We present the Einstein frame potential in the i?^-corrected f{R) model and consider 
its characteristics analytically. In Sec. Ill, we derive the analytic solutions of the inflation and reheating in the 
Einstein frame. We shall adopt the slow roll and the fast roll approximations and solve the field equations in each 
era. Sec. IV contains the results of the numerical calculation. We confirm that the analytic solutions are sufficiently 
in agreement with the numerical results. We also consider the behavior at the end of the reheating. In Sec. V, we 
consider the connection between observables and model parameters and discuss its allowed ranges. Sec. VI is devoted 
to conclusions and discussion. Throughout the paper, we adopt units c = h = 1. 



II. BASIC EQUATIONS 

We start to review the basic equations of f{R) gravity in both the Jordan frame and the Einstein frame. To avoid 
confusion of the frames, we fix the subscript J and E to physical quantities in the Jordan frame and the Einstein 
frame, respectively. Otherwise we declare the frame in which the quantity is defined, /(i?) gravity is defined by the 
action 



S ^ d x^-gj 



f{Rj) + CM{9') 



(1) 



where Cm is the Lagrangian density for the matter sector and M-p\ is the reduced Planck mass. By varying the action, 
we obtain the field equation in the Jordan frame 



Ri,F{Rj) - \giJ{Rj) + {gi^U - V^V,)F{Rj) = ^ 



with 



F{Rj 



_ df{Rj 



dRj 



(2) 



(3) 



We regard the Jordan frame as the physical frame. However, for our purpose analyzing the infiation and the 
reheating in f{R) gravity, the formulation in the Einstein frame is useful because it contains the additional degree 
of freedom more explicitly as a scalar field and enables us to use the analogy of the single-field inflation. We can 
recast the theory to the Einstein gravity with a scalar field by choosing the conformal transformation of the metric 
as gj^^ = F{Rj)g^j^. The canonical scalar field is defined by 



FiRj) = e 

By the conformal transformation, the action is rewritten as 



S 



with the potential term 



PI 



Ri 



,(j)-V{<j,)+CM\e 



V{4>) 



Ml, Rj{4>)F{Rjm- f[Rjm 



2 F{Rj{4>)f 
Then, the Einstein equation in the Einstein frame reduces to 



Hi 



dHE 
dtK 



1 



2A/2j 



dtt 

d<j) 
dtE 



o 17^ +V{<I^)+PE 



Pe 



Pk 



(4) 
(5) 

(6) 

(7) 
(8) 
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The equation of motion for the scalar field is 

From the conformal transformation, the time and scale factor in both frames are connected by 

1 <^ 1 

dtj = e -^^'pidtE, aj = e v^^pia^. (10) 

The transformation of the Hubble parameter is derived from the above definitions. 

1 <*> f 1 d6\ 
Hj = ev^^^pi He 7= r— . 11) 

By definition in Eq. (3), the energy momentum tensors of the matter sector in both frames are connected by 



^ e y'^~iT-^ . (12) 
For perfect fluid Tj^ — diag(— p, P, P, P) in each frame, the energy density and the pressure are related as 

P£ = e'^^WTpj, Pe ^ e^^^T^iPj. (13) 

Note that the energy density in the Einstein frame couples with the scalaron. 

In the inflation and reheating in /(i?) gravity, there is no inflaton field from the point of view in the Jordan 
frame. Consequently, particle creation occurs not through the decay of the inflaton but through the gravitational 
reheating [4-8]. Let us consider the gravitational particle creation in the Jordan frame. We introduce a minimally or 
nonminimally coupled massless scalar field Xj which describes the created particles, into the matter action, 



S = d xy/-gj 



„,„ . 1 „ 1 



(14) 



Since the scalar field x couples with the metric in the Jordan frame, the radiation (massless scalar particle) is created 
purely gravitationally. Adopting the standard treatment of the quantum field theory in curved spacetime, we can 
expand x in Fourier modes with the annihilation and creation operators. Then, computing the Bogolubov coefficients 
in the expanding Universe, we obtain the number density of the created scalar particles [4-7] , 

The above equation holds regardless of the functional form of f{R). Note that the particle creation is sourced by Ricci 
curvature. Since the Ricci curvature is significantly suppressed during the reheating era as we shall see below, particle 
creation hardly occurs. On the other hand, inflationary dynamics is the same as that of the inflation. Therefore, 
we can use the approximated formula for the model and turn off the particle creation during the reheating era 
when we perform numerical calculation. In the R^ model with f{R,j) = Rj + i?j/(6M^), the energy density of the 
created particles is 



(10) 

where denotes the relativistic degree of freedom. In the present paper, we consider minimally coupled massless 
scalar field and set ^ = hereafter. The evolution equation for the energy density of radiation is then 

|^^-4i/.p, + ^. (17) 
at J ll527r 

The pressure is obtained from the energy conservation equation, 

p^^fU_jMR^_ (18) 
3 'SAdGnHj ^ ' 
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Finally, we introduce gR^-AB model [34] 
present Universes, 



which describes the accelerated expansions in both the early and the 



f{Rj) = {l-g)Rj+gMH\og 



cosh{Rj/SAP - b) 
coshfe 



6Af2' 



(19) 



where g, &, 5, and M are positively-defined model parameters. 5 describes the ratio between the energy scale of the 
present DE to the primordial DE and takes dramatically small value. M will be fixed to M/Mp\ 1.2 x 10~^ later 
in Sec. V by the temperature fluctuations of the cosmic microwave background (CMB) anisotropy. g is constrained 
in the range of < 5 < 1/2 by the stability conditions of f{R) gravity: F{Rj) > and dF{Rj)/dRj > 0. Moreover. 
g and b arc further constrained by the other stability condition as we shall show in Sec. V. The function F{Rj) is 
given by 



FiRj) = l-g 



Rj 
MP 



gta.iih{Rj/Ap6 -b). 



(20) 



The above ,f{R.j) function is equivalently written in the following form: 



f{R:> 



Rj 



gM^S log 



1 



i?vac = 2gA'P6 {b + log(2 cosh b)} . 



3_ 

6M2 



(21) 
(22) 



In Eq. (21), the fourth term dominates at high curvature regime Rj > and causes inflation [3]. The third term 
alters the reheating dynamics after the inflation, which is characteristic of the gR^- AB model. The second term plays 
the same role as the current cosmological constant. By substituting the action into Eq. (2), the equation of motion 
in the Jordan frame is given by 



12 



^—IaPh]-^m\h'j 

2 2 



cosh{Rj/M'^5-b) 



cosh 6 



H"]) tanh 



Rj 
APS 



- b 



3g{HjHj + 4H'jHj) 
6 cosh^ {Rj /APS - b) 6A/^i ' 



AP 

— tpj 



(23) 



The prime denote time derivative with respect to tj. As expected from the specific form of f{Rj) function, the 
equation of motion in the Jordan frame is nonlinear differential equation and quite complicated to solve. 

The dynamics of inflation and reheating can be more intuitively understood from the potential in the Einstein 
frame, which is depicted in Fig. 1 (a). From the deflnition of Eq. (6), we can interpret the shape of the potential 
in the following way. First, we notice that from Eq. (20), F{Rj) is almost a step function at Rj/AI^ ~ bS with the 
change of the amplitude from = 1 — 2g to 1. In terms of the scalaron, the step corresponds from (f)/AIp\ = \/61og7 
to 0, where 7 = \/l — 2g. In other words, during the scalaron moving in the range V^log7 < 0/Mpi < 0, Rj 
remains almost constant value, Rj/AP ~ bS. Outside this interval, F{Rj) is approximated by F ~ 1 + R/iA'P for 
Rj/hP > bS, I.e., 0/Mpi > 0, and F ~ 1 - 2g + R/SM"^ for Rj/AP < bS, i.e., (^/Mpi < VGhg-f, respectively. By 
using these approximations, we can derive the potential analytically in terms of (/): 



Aq^Ap 



Ye V 3 Mpj 



(0/A/pi > 0) 
(</./A/pi < V61og7) 



(24) 



In these two regions, the potential is the same as that in the R^ inflation. On the other hand, the characteristic 
plateau shows up for -\/61og7 < (p/AIpi < 0. By using Rj/AP ~ bS, we can approximate the potential as 



V{4>) 

M'^yAP 



bSe^ 



7?r - f{bSAP)/AP 



(V61og7 < 0/Mpi < 0) 



(25) 



2e/ 



From this expression, we can estimate the height of the bump in the plateau. We can obtain the position of the local 
maximum by solving V' — 0. The solution is (^/Mpi ~ ^3/2 log[2(l — 2g) + bS/3\ = 4>m/Mp\. The potential has the 
local maximum when (j)^ satisfies •\/61og7 < ^^A^pi < 0, namely, (3 + 6i5)/12 ^ 5 ^ (3 + 6(5)/6. For 5 <C 1, this 
condition amounts to 1/4 < g < 1/2. Therefore, in Fig. 1 (c), the potential for g = 0.35 possesses the local maximum 
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Figure 1: Inflaton potential of the gR^-AB model in the Einstein frame, (a): Potential for g = 0.35, b = 5, 5 — 5 x 10~^. For 
(f) > and (fi < %/6AfpiIog7, the potential is similar to that in pure inflation. On the other hand, for •\/6Ai/"pi log 7 < (f> < 0, 
there is the characteristic plateau. The scalaron starts slow rolling from (j) > 0, and enters the plateau with large kinetic energy 
and oscillates inside it. (b) - (d): How plateau changes its shape when one parameter is changed from the value in (a). The 
typical height of the plateau is bgS from the analytic estimation. 



and the false vacuum. However, as g approaches g = 1/4, the local maximum is likely to disappear. The potential 
heights at the right edge, the local maximum, and the left edge are given by 



bgS, 



(1-25)' 



F(V6Mpilog7) 



0, 



(26) 



where we used 6 <ti I and log(cosh6) ~ b. These estimations well explain the parameter dependences of the potential 
shape in Fig. 1 (b) - (d). 

Next, let us consider the evolution of the scalaron. The scalaron starts slow rolling from (f) > and plays a role of 
the inflaton. For > 0, the potential is the same as that of a pure model and is almost independent of model 
parameters b and S. Thus, the scale factor in the Einstein frame experiences quasi-de-Sitter expansion. In this 
case, the scale factor in the Jordan frame also evolves exponentially, because the amplitude of (j) slowly varies and 
thus the scale factors in both frames are related by multiplying an almost constant factor. As the scalaron approaches 
(/> = 0, it rolls faster and enters the potential plateau with the kinetic energy larger than the potential energy. Then 
the scalaron oscillates in the plateau, gradually loses its kinetic energy, and finally reaches the false vacuum at = 
because the chameleon effect lifts up the potential when the energy density of matter does not negligible [35] . During 
this oscillation, the scale factor in the Jordan frame undergoes the periodic evolution due to the exponential factor in 
Eq. (10). We shall see these situations in the next section. 
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III. ANALYTIC SOLUTIONS 



To investigate the reheating dynamics in the Jordan frame, we work in the Einstein fr-ame and derive the analytical 
expressions for the motion of the scalaron in the inflaton potential. For (j> > 0, the potential is almost reduced to 
that of a pure model, then we can use the slow-roll approximation in Sec. Ill A and derive the analytic solutions. 
After the end of the inflation, the kinetic energy of the scalaron is dominant. Therefore, we can neglect the small 
structure of the potential during the oscillation in the plateau and use the approximation, Rj ~ bS. We shall explore 
this case in Sec. IIIB. Since in both the slow-roll and oscillation regimes the energy density of the created radiation 
is subdominant in comparison with the energy density of the inflaton, we neglect its backreaction to the background 
dynamics to derive the analytic solutions. However, the energy density of radiation eventually becomes the same 
order as the total energy of the inflaton, and the reheating ends. In the following, we define dimensionless variables 
i = Mt, <j) = 0/A/pi, V = V/APMly, H = H/M, R = R/AP, p, = p^/M^iAP and abbreviate the hat in this section 
to avoid the complexity of notation. 



A. Slow roll approximation 



The inflaton starts to roll down the potential at = tEMi {tj = ijani) with small kinetic energy compared to the 
height of the potential. The energy density of the radiation is also negligible. The field equations (7) - (9) with slow 
roll approximations are 



He 
He 



Vi4>) 

3 ' 



(27) 
(28) 
(29) 



where the dot denotes the derivative with respect to the time in Einstein frame and the potential is approximated by 
Eq. (24). We set the initial condition as = (f>ini- Other initial quantities are fixed from the above equations with the 
approximated potential. By substituting the potential into the field equations, we can derive the following solutions: 



HE{tE) 

aEits) 

Rj{tE) 

where we use dini = e ^ instead of (b 



dt 



-{tE — tE.m 



log 

2 



: 3(e\/f'^(*'=) 
itself. 



itE-tE,ini)/2 
1), 



3/4 



(30) 
(31) 
(32) 

(33) 
(34) 



The remaining task is to write down in terms of tj and to convert the above solutions. Jordan frame time is 
derived by integrating Eq. (10), 



tj{tE) = tjM 



i\ldl,,^ - {tE-tEAr. 



3di, 



and its inverse function is 



tEitj) — tE,ini + {tj — tj,ini) 



rfini — g(t,7 — ijani) 



(35) 



(36) 



Thus, substituting Eq. (36) into the solutions (30) - (34) and using Eqs. (10) and (11), we obtain the analytic solutions 
in terms of the Jordan frame quantities. 
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Next, we derive an analytic solution for pj by performing integration in Eq. (16) with Eqs. (10), (30), (33), (34) 
and (36), 



Pr{t,j 



9* ( M V 1 



11527r \MpJ 165*2 

-2a/35(454 - 1452 + 15) + 3e'S' |-6di„ie-3''i"i (12<i - 14<i + 5) + 5^3^ (erf(V3di„i) + erf (5)) } 



(37) 

where S = {tj — tj^i^i — 3(iini)/V3 and erf(x) is the error function. Hereafter, we refer pj as pr, especially denoting 
radiation component. 

Finally, we focus on boundary conditions. By using the above analytic solutions, we define the time Ie = tEo or 
tj = tjQ when the inflaton reaches = for the first time. Strictly speaking, when ~ 0, the slow-roll approximation 
is not valid anymore and the boundary conditions are very sensitive due to the the sudden transition of the potential 
at (/) ~ 0. Therefore, we should use the boundary conditions obtained from numerical computation. We only rely on 
the analytical boundary values as the estimator. We shall revisit this point in Sec. IV. 

The time Ieo is analytically estimated by using the analytic solution (30), 

3 

tEO = tE,ini + 2^*^™' ~ '-"^^-^ 



In terms of Jordan frame time. 



The boundary conditions are 



From Eq. (7), 



tjO=tj,ini + 3(di„i-l). (39) 



ho = 0, 00 = -^! ai5o = ais,i„irfi„?e^('^'"'-il (40) 



Heo = I, Hjo = I (41) 



The energy density Pr at tj = tjo is given by setting S = — \/3 in Eq. (37) and taking large (pini limit, i.e., dini ^ 1, 



Pro = 



co9* f A/ ^ ^ 



11527r V A/pi 



The coefficient cq is analytically given by [l8 + 5e'^-\/37r (l — erf(-\/3))] /16 w 1.4, but at this intermediate stage from 
the slow roll to the fast roll, the slow-roll approximation does not give the precise value. So we will use the result of 
the numerical calculation for the precise determination of cq later. We notice that since p^o ^ jq for typical values 
5* = 106.75 and M <C Mp\, we can neglect the energy density of the created radiation for the purpose of deriving 
analytic formulas of reheating dynamics in the next subsection. 



B. Fast roll approximation 

When the inflaton falls down to the plateau on the bottom of the potential that lies a/6 log 7 < < 0, its kinetic 
energy is much greater than the potential energy, which is typically of the order of bgS. Moreover, it is also greater 
than the energy density of radiation. Thus, we will use two approximations to derive the analytic solutions. First, 
the energy density created at the slow roll regime is not dominant as wc mentioned above. Second, particle creation 

we can use the fast-roll approximation in Eqs. (7) - (9): 



is negligible during the fast- roll regime because the source term i?j keeps an almost constant value ^ 1. Thus, 



(42) 



a2 



He = -y, (43) 
+ 3i?£;0 = O. (44) 
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In the plateau, the inflaton repeatedly oscillates between the left and right walls of the potential. We separately 
consider the time intervals dependent on the direction of the motion of the inflaton and derive the analytic solution 
in each regime. 

The inflaton is reflected at the left side of the plateau at = \/6\ogj. We regard the reflection occurs instantly 
and define the first refiection at t = tEi- After that, the inflaton reaches the right wall at = and is reflected at 
t = tE2- Thus, we can periodically define tEn until the inflaton stops somewhere. 

First, let us consider the interval Ieo < tE < tEi- Eliminating from Eq. (42) and (43) and then solving it, we 
obtain the Hubble parameter and the scale factor 

aEits) = ai^o[3i^Bo(^i^-^Bo) + l]'^^ (46) 

From Eq. (42), we obtain (f> ~ ±\/6He- For t.Eo < Ie < tEi, since the inflaton moves in the left direction, we choose 
(j) = —\/6He- Therefore, the solution is given by 

HtE) = -^\og [SHEoitE - ifio) + 1] • (47) 
Next, we move to the interval tEi < tE < tE2- By using cp — +\^He, we can derive 

aE{tE) = aEi[iHEiitE-tEi) + l]^^\ (49) 
cpits) = ^J^\og[SHEi{tE-tEi) + l] + V6logj . (50) 

The matching conditions between the two intervals are determined by setting (f>{tEi) = ■\/61og7 in the solution for 
tso <tE < tEi as 

tEl = *B0 + ^T7 ' (51) 

Hei = HeoI^ (52) 
a-Ei = aE07^^- (53) 

By substituting them into the solutions (48) - (50), we obtain the same expressions as in Eqs. (45) and (46) for 
-H^s(^s) and asitE)- However, 4>{tE) is different from Eq. (47). It is given by 

^{tE) = y|log [iHEoitE ~tEo) + l\ + 2\/61og7. (54) 
Likewise, by solving the recurrence equations, we conclude that the boundary conditions for general n are 

(55) 

(56) 
(57) 

(58) 

Especially, it is noteworthy that the following relation holds regardless of the parity of n: 

SHEnitE - tEn) + 1 = t'^""'^ [3i?£o(ti=; ^ tfio) + 1] • (59) 



tEn 






3Heo 




HEn 






^En 


= aso7~": 








(n : even) 


<t>En 


\ \/61og7 


{n : odd) 
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Thus, the solutions for tEn-i < ts < tEn are 

Heo 



HE{tE) = 
aE{tE) = 



'iHEn{tE — tEo) + 1 ' 

aEo [3i/Bo(i£; - teo) + 1]^^^ , 

-y|log [SHEaitE - tEo) + 1] - (n - l)V61og' 
log [3HEQ{tE - i_Eo) + 1] + 7i\/61og7 



{n : odd) 
[n : even) 



{n : odd) 
(n : even) 



-V6HEitE) 
VEHsitE) 



(60) 
(61) 

(62) 
(63) 



Notice that the solutions for (j) and (f) have different expressions for the different parity of n because of its direction of 
the motion. 

Once the solutions in the Einstein frame are at hand, it is straightforward to convert them to those in the Jordan 
frame. The Jordan frame time tj for tEn-i < tE < tEn is obtained by integrating e~'^/^. 



Tl-l r 



tj{tE) 

where tjn is given by 




{3HEo{tE — tEo) 
{3HEo{tE — tEo) 



,4/3 



-4(n-l) 



1} -' ^ - 7 

;L 1.2/3 _ ^-2(n-l) 



Jn 



JO 



tjo 



(74 + 72 + 2)(7-3(«-i)-l) 7-'^-l 



4Heo{j^ + 1^ + I) 



(7^ + 7^ + 2) (7 



-3 11 



1 



{n : odd) 
(n : even) 

(n : odd) 
(n : even) 



(64) 



(65) 



4-ffiso(7^ + 7' + l) 

From Eqs. (10), (11), (60), (61), and (64), the Hubble parameter and the scale factor in the Jordan frame evolve as 

273("-i)frB0 



Hjit, 





^ J ajo7 



473("-i)iJso(ij-tj. 



1 



{n : odd) 
[n : even) 



-(n-l) 



4-7 



3(ri-l) 



HEo{tj — tjn-l) + 1 



1/2 



(n : odd) 
(n : even) 



(66) 



(67) 
(68) 



Hj(tj) periodically takes discrete behaviors at tj = tj„ and vanishes for even n. However, to be precise, Hj is not 
exactly zero because we here neglect the contribution from the energy density of created radiation and the potential. 
We shall numerically confirm this point in the next section. We define Hjn by lim^^±Q H(tjn + e) where ± for even 
and odd n, respectively. From the Hubble parameter (66), it is given by 



Hjn — 2He,i — 7'^"iJjo- 



(69) 



Here we used Eqs. (41) and (56). 

Let us remark the time averaged behavior during the oscillation. Since Eq. (55) implies \og(tEn — tEo) oc n, the time 
intervals when inflaton goes to the left and the right are equal in terms of the logarithmic Einstein-frame time. On 
the other hand, Eq. (64) yields log(tj — tjo) — p\og{tE — tEo) with p = 4/3 and 2/3 for the left-directed regime and 
the right-directed regime, respectively. Hence, the left-directed regime is twice as long as the right-directed regime in 
terms of the logarithmic Jordan-frame time. Taking care of these facts, we can derive the average power of the Jordan 
frame quantities. For instance, as for time duration, since the average power is (4/3 + 2/3)/(l + 1) = 1, we can say 
tE is proportional to tj in average, i.e., {tj — tjo) ex. ts — tEo- Since logoj ~ qlogtj with q = 1/2 and for the 
left-directed regime and the right-directed regime respectively, the averaged power is (l/2x 2 + Ox l)/(2 + l) = 1/3, 
namely, {aj{tj)) oc {tj — tjo)^^"^. Hence, the averaged Hubble parameter is written as 



(Hjitj)) 



H 



JO 



JO) 



(70) 
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Integrating this gives the averaged scale factor 



Mtj)) 



a. JO [iHjo{tj-tjo) + l] 



1/3 



(71) 



In summary, neglecting the small structure of the plateau and the energy density of the radiation enables us to 
solve the differential equations analytically. However, in the final phase of the reheating, they become important. We 
shall consider this effect in the next section. 



As we derived in Sec. II, the typical height of the local maximum of the potential plateau is Vmax d. The oscillation 
regime ends when the kinetic energy of scalaron is of the same order of Vinax- Thus, 6 determines the end time of the 
reheating. However, since S is the ratio between the energy scales of DE and of inflation, it has to be dramatically 
tiny value ^ 10~^^°. Therefore, it is difficult to carry out numerical calculation for realistic model parameters. In this 
section, using the modest value of 6, we confirm the validity of the analytic solutions obtained in the previous section 
by comparing them with numerical results. Then we extrapolate the analytic solutions to the end of the reheating 
and estimate the reheating temperature as a function of the model parameters. Note that hereafter we explicitly fix 
the hat to normalized dimensionless physical variables, which was abbreviated in the previous section. 



We performed numerical calculation in the Einstein frame, i.e., we solved the coupled evolution equations: Eqs. (8), 
(9), and (17). In particular, to avoid the complexity due to nonminimal couplings in radiation and matter sectors, we 
use Eq. (17), instead of the continuous equation in the Einstein frame, by translating it to the Einstein frame with 
Eqs. (4), (10), and (11). After the numerical calculation, we obtain physical quantities in the Jordan frame by the 
inverse conformal transformations. 

Figure 2 illustrates the numerical results and the analytic solutions for a model parameter set: g = 0.35, b = 
5, (5 = 5 X 10~^, M/Mpi = 1.2 x 10^^. Since observationally required value for 6 is too tiny to perform the numerical 
calculation, we chose this S just for demonstration and compared its results with the analytic solutions. As mentioned 
before, the slow-roll approximation does not hold at the intermediate stage from the slow roll to the fast roll. So we used 
the boundary conditions obtained by the numerical calculation as the initial conditions for the fast-roll approximated 
solutions: cq = 0.72 and Hjo = 0.26 (These are different from the values obtained from slow-roll analytic solutions 
by about a factor of 2). The analytic solutions well approximate the numerical results. In the slow-roll regime, the 
scale factor undergoes quasi-de Sitter expansion. In the fast-roll regime, the inflaton oscillates two times inside the 
potential plateau. We see that the scale factor evolves as a,j{t,j) oc (tj — tjoY^'^ and aj(tj) ~ const in order. Ricci 
scalar stays constant value ~ bS except for several spikes during the reheating. In contrast to the analytical solution, 
Hj does not vanish when the infiaton moves to the right direction in the potential plateau because of the contribution 
from the energy density pr of the created radiation. At late time, the fast-roll approximation becomes worse as the 
inflaton loses its kinetic energy. Finally, the infiaton reaches a false vacuum at = 0. 



Since the above calculation have performed for (5 = 5 x 10~*, the scalaron ends the oscillation quickly before 
radiation dominates the Universe. However, the scenario is different for S ^ 10" ^^'^: the reheating appropriately ends 
by radiation domination. This is because V{0) and Vinax is too small to compete with the inflaton kinetic energy. For 
the following, by using the analytic solutions, we estimate both the times when the scalaron stops the oscillation and 
when radiation dominates the Universe. 

First, let us estimate when the oscillation stops. It is estimated from the equality time of the kinetic energy of the 
scalaron and the local maximum of the plateau: (fP/2 ^ Vmax- We use the analytic solution for cj) in Eq. (63) and the 
asymptotic behavior (He) ^ l/3(t_E ~ i_Eo) ~ l/3(t,/ — tjo). Thus, the scalaron ceases the oscillation at 



where we used Knax by Eq. (26). Substituting 6 lO-^^o^ ^.j ^ 2. 2 x IQ-'^Mpi, and 5,6 - 0(1) yields tjs ~ 
10^^ GeV~^, which is close to the Hubble time Hq^. To be precise, we should take the radiation and matter 



IV. NUMERICAL RESULTS 



A. Comparison with analytic solutions 



B. End of the reheating 
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Figure 2: Evolution of inflaton, e-folding number, Hubble parameter and Ricci scalar for model parameter g = 0.35, b = 5,S — 
5 X , M/Mpi = 1.2 X 10"^. Numerical results (blue, solid), slow roll analytic solution (magenta, dashed) and fast roll 
analytic solution (green, dot-dashed) are presented. 



dominated epochs into account in the following evolution of the Universe, but it does not change the conclusion that 
the scalaron oscillation continues for the order of the Hubble time. 

Next, let us estimate when the radiation dominates the Universe and the reheating ends. We compare the energy 
density of radiation due to particle creation and that of gravity. Since Rj ~ bS during the oscillation in the plateau 
of the potential, we can expand F around Rj = bS and take its linear order only. 



+ {Rj-bS)^eVii. 

Thus, we can explicitly connect the Ricci scalar in the Jordan frame with the inflaton as 

35(e\/i<^(*^) -l + g(6+l)) 



Rjiti 



3.9 + (5 



(73) 



(74) 



As we mentioned at the beginning of Sec. IIIB, the particle creation during the plateau oscillation phase is negligible 
because Rj ~ 6(5 ^ 1. Therefore, pr is approximately given by 



{Pr{tj)) = Pro 



a, J a 



Pro 



(75) 



where we used Eq. (71). As for gravitational contribution, it is convenient to define the effective energy density of 
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Figure 3: Evolution of energy density of radiation pr (blue) and effective energy density of gravity pg ~ 3Hj (red) for model 
parameter g — 0.35, 6 = 5, 5 = 5 x 10~*, M/Mp\ = 1.2 x 10~^. For energy density of radiation, the analytic solutions for slow 
roll regime (magenta, dashed) and fast roll regime (green, dot-dashed) are also presented. 



gravity by the equation of motion in the Jordan frame in Eq. (23): 
1 



Pg 



{Pt+ Pg) 



3M, 



M2 



' gM'^H] - 2H"jHj + {h'jY - 6H'jH^j + gM^{H'j + Hj) tanh 
cosh(i?j /M^6 - b) ] GgiHjHj + AHjHj) \ 



Rj 
M^6 



6 



coshfe 



When Pj. is negligible compared to Pg, from Eq. (70), the energy density of gravity is reduced to 

1 



{Pg 



m.j? 



(76) 



(77) 



(78) 



Figure 3 represents the evolution of pr and pg « 3i/j. Throughout the inflation and reheating, the energy density 
of the radiation is subdominant for our choice of (5 = 5 x 10^^. The radiation is mainly created at the end of inflation 
and decays as Pr oc a~^. Since the scale factor evolves as t^/^ and constant periodically, Pr correspondingly evolves 
as tj^ and constant during the left-directed and right-directed regimes, respectively. As wc mentioned, 3Hj does 
not vanish because of the contribution from the radiation. At Mtj ~ 50 and 1000, 3Hj is smaller than p^. This is 
because the required step width at the abrupt transition is so tiny that Hj over-decreases. This tiny discrepancy is 
unimportant and docs not affect the subsequent evolution. 

The end time of the reheating and the reheating temperature Ti-ch is defined by the condition (pr) = (pg)- Using 
Eqs. (75) and (78), wc obtain 



^ J,reh 



-JO 



JO 



.3/2 
PrO 



3.8 X 10^(co.9*)"^/^H 



/ M 



PI 



and 




« 9.6 X 10-3 (co5. 



Nl/2 



M 



(79) 



(80) 



Since these equations include Hjq and pro, the end of the reheating is sensitive to the boundary conditions at the 
transition from the slow-roll regime to the fast-roll regime. For = 106.75, M/Mp\ = 1.2 x 10^^, and numerically 
determined values, Cq = 0.72, Hjq sa 0.26, the time and temperature of the reheating arc tj^rch ^ ^ jo ~ 2.2 x 10^^ and 
Treh ~ 3.0 x 10^ GeV, respectively. 
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V. CONSTRAINTS ON THE MODEL PARAMETERS 



111 this section, based on our analytic solutions and numerical results obtained in the previous sections, we discuss 
the allowed ranges of the model parameters. First, we present the constraint on the energy scale M in Sec. VA. 
It is fixed from the amplitude normalization of the CMB. Second, we consider the other parameters, b, and S in 
Sec. VB. Once M is determined, we can predict the other observable of the inflation and discuss its consistency with 
observations, g, h, and 5, cannot be constrained from observational data because the reheating temperature in Eq. (80) 
is independent of those parameters. However, since the gR^-AYi model must realize current cosmic acceleration, its 
magnitude and stability constrain the allowed range of g, h, and 5. 



A. Constraint on M 



When Rj ^ A-P , the gR^-KQ model (21) can be approximated to R^ inflation, in which the primordial spectrum 
of a scalar mode at the leading order in the slow-roll parameter ei is given by [36] 



Vs 



1 



M 



(81) 



where ei = —Hj/Hj. This slow- roll parameter is related to the e- folding number between the end of inflation and the 
horizon crossing of the mode whose comoving wave number k corresponds to the CMB scale today. From the analytic 
solutions during the slow-roll regime, Eqs. (11), (30) - (32), the Hubble parameter in the Jordan frame is written as 



Hjitj) 



3t(£j)-1 



6Ti/2(f^) ' —'-"mi 3 

For dini ^ 1, expanding Eq. (82) in powers of tj around ijjni gives 

1 



^ini ^ (^ J ^ J,ini ; 

6 



HjSni ~ 7(^.7 ^ ij.ini) , 
6 



(82) 



(83) 



where Hjini = (iini/2. The e-folding number between the end of inflation at tj,cnd and the horizon crossing of the 
CMB mode at tj^ is given by 



*.7,cnd 



Hidt 



1 



2ei(t 



Jk) 



(84) 



where Hjk = Hj{tk) and we used the fact that Hjj^ is constant when we performed the integration. Then Eq. (81) 
is expressed as 

2 



Vs 



Nl ( M 



(85) 



247r2 \My>\ 

Using Eq. (75), we obtain the e-folding number when the comoving scale of CMB crosses the horizon during inflation 



TVfe « 66.2 - i log 



1-a. 

0.7 



-ilog 



Co \ / a* 



0.72/ V 106.75 



(86) 



Since the parameters il„i and hardly change Nk, we set it to Nk = 66. From the temperature fluctuation of CMB 
anisotropy [9], the amplitude of the power spectrum, Vs = (2.445 ± 0.096) x 10~^ at = 0.002 Mpc^"^, fixes the 
parameter M to 

M 

— « 1.2 X 10"^ . (87) 
Mpi 

At the CMB scale, the spectral indices of the scalar and tensor modes and the tensor-to-scalar ratio are given by [36] 



ns - I 



dlogVsik) 



d\og k 
dlogVrik) 



d\og k 
Vt 12_ 



k=aH 



k=aH 



2 

'^k ' 

3 



(89) 
(90) 
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at the leading order in the slow-roll parameter. For the above choice of Nk, ns ~ 0.97 and r k, 2.8 x 10 which are 
consistent with observational bounds [9] . 



B. Constraints on g, b, and S 



The parameters g, b and S considerably alter the dynamics of the reheating in the gR^-AB model. However, as seen 
from Eqs. (75) and (80), the reheating temperature and the radiation energy density at that time does not depends 
on these parameters. So the constraint on g, b, and S comes not from the CMB observation but from a stability 
condition of a de-Sitter vacuum. 

From Eq. (2), 



3nF{Rj) + RjF{Rj) - 2f{Rj) ^ SnGTj , 



(91) 



where Tj = g'^j^T^^ and □ = (l/^/^gj) [\/~9J For the existence of a stable solution of a de-Sitter vacuum 
{Rj = const., Tj = 0), the following equation has to be satisfied: 

RjF{Rj) - 2f{Rj) = . (92) 

For Rj ^ AP, substituting Eq. (21) into Eq. (92) and using i?vac ~ igbAPS lead to the equation 



Q{y) = y--igb + 2g 



log 1 + e 



-2{y-b) 



1 + e^^y-^^ 



. 



(93) 



where y = Rj/M'^6. This function Q{y) typically has the shape shown in Fig. 4. Therefore, a stable de-Sitter vacuum 
exists if Q (yo) = has the solution y = yo > 1 at which Q (yo) > and Q{yo) < 0. The boundary of the allowed 
parameter region of b and g can be obtained by solving Q{yo) ~ and Q (yo) = under the condition Q (yo) > 0. 
We cannot solve the above equations analytically. Instead, we fit the numerical solution and obtain the allowed region 
for g as 



1 0.28 1 

4 ^ (6-0.46)0-81 - -"^ - 2 



(94) 



This region is shown in Fig. 5. 

Once the parameters AI, g, and b are fixed, S should be determined so that the current observation of accelerating 
expansion is reproduced. With the Ricci curvature of the present universe, i?vac ^ 10^i^"i\/pi, the parameter S is 
given by 



1 i?v 



2gAP{b + log[2 cosh 6]) Agb AP 



(95) 
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Figure 5: Allowed parameter region of g and b. Numerical solution (solid curve), fitting (dotted curve). Above this curve, 
stable de-Sitter solutions exist. 



VI. CONCLUSIONS AND DISCUSSION 

We have studied the inflation and reheating dynamics in f{R) gravity, especially in gR^- AB model. This model is 
capable to describe both accelerated expansions in the early Universe and the present time. In the Einstein frame, 
the inflaton potential of this model possesses a plateau and a false vacuum in the bottom of the potential. These are 
different features from original inflation model and they significantly change the reheating process. We have derived 
the analytic solutions in the slow-roll inflation regime and the fast-roll oscillation (reheating) regime. We have also 
carried out the numerical computation including the backreaction from particle creation, and have confirmed that both 
results agree well. According to the existence of the potential plateau, the particle creation via gravitational reheating 
mainly occurs in the slow-roll regime and is inefficient during the the fast-roll oscillation regime. Consequently, in 
contrast to the R^ inflationary scenario, the reheating era lasts longer. Another interesting feature of this model is 

1/3 

that the averaged time evolution of a scale factor is proportional to tj because of the periodic abrupt changes of the 
Hubble parameter. 

Based on these results obtained from our analytic and numerical calculations, we have given the constraints on 
the model parameters. The parameter Af is pinned down by the observational amplitude of CMB temperature 
fluctuations. Also the value of 5 is selected to correctly reproduce the current accelerated expansion of the Universe. 
On the other hand, the parameters g and b arc poorly constrained because these parameters affect only the reheating 
dynamics after the inflation. To more tightly constrain g and 6, we need observations that can probe at much smaller 
scales than those of CMB and large-scale galaxy surveys. In the future searches for primordial black holes and 
the direct detection experiments of gravitational waves would provide new observational windows for the reheating 
dynamics in modified gravity. 
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